import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import minimize

def calculate_y1(x1, x2, x3):
    y_i = 0.496 * x1 + 0.567 * x2 - 100.046 * x3 + 5.499 * x3**2 - 0.090 * x3**3 + 1644.941
    return y_i
def calculate_y2(x1, x2, x3):
    y_i = -0.285 * x1 + 0.053 * x2 - 3.721 * x3 + 0.224 * x3**2 - 0.004 * x3**3 + 105.786
    return y_i
def calculate_y3(x1, x2, x3):
    y_i = -4.320 * x1 + 0.362 * x2 - 10.379 * x3 - 0.101 * x1 * x3 + 0.047 * x1**2 + 0.395 * x3**2 - 0.001 * x3**3 + 189.837
    return y_i
def calculate_y4(x1, x2, x3):
    y_i = 22.863 * x1 - 1.597 * x2 + 32.268 * x3 - 0.172 * x1 * x3 - 0.401 * x1**2 - 2.368 * x3**2 + 0.044 * x3**3 + 73.167
    return y_i
def calculate_y5(x1, x2, x3):
    y_i = -1214.089 * x1 - 182.990 * x2 - 69.051 * x3 + 0.494 * x1 * x2 + 2.878 * x1 * x3 + 0.494 * x2 * x3 + 42.671 * x1**2 - 0.585 * x3**3 + 0.723 * x2**2 + 0.363 * x3**2 + 23669.459
    return y_i
def calculate_y6(x1, x2, x3):
    y_i = -0.048 * x1 + 0.006 * x2 + 0.068 * x3 + 1.442
    return y_i
def calculate_y7(x1, x2, x3):
    y_i = -0.012 * x1 + 0.108 * x2 + 0.493 * x3 + 140.43
    return y_i
def Hot_wet_comfort_performance(m4, m5):
    return -((m4-183.5229167)/73.27417756+(m5-2626.767 )/402.387)
# 定义目标函数
def f1(x):
    x1=x[0]
    x2=x[1]
    x3=x[2]
    return Hot_wet_comfort_performance(calculate_y4(x1, x2, x3),calculate_y5(x1, x2, x3))


def f2(x):
    x1 = x[0]
    x2 = x[1]
    x3 = x[2]
    y1=calculate_y1(x1, x2, x3)
    y2=calculate_y2(x1, x2, x3)
    y3 = calculate_y3(x1, x2, x3)
    y4 = calculate_y4(x1, x2, x3)
    y5 = calculate_y5(x1, x2, x3)
    y6 = calculate_y6(x1, x2, x3)
    y7 = calculate_y7(x1, x2, x3)
    return (y1-1331.605625)/233.5740652+(y2-94.32432083)/8.161031342+(y3-117.27)/26.17596577+\
           (y4-183.5229167)/73.27417756+(y5-2626.767 )/402.387+(y6-2.07)/0.84 +(y7-160.03 )/6.63


# 定义ε约束
def epsilon_constraint(x, epsilon):
    return epsilon - f2(x)


# 定义多目标优化的主函数
def epsilon_constraint_method(f1, f2, bounds, epsilons):
    solutions = []
    for epsilon in epsilons:
        # 约束条件
        cons = [{'type': 'ineq', 'fun': lambda x, eps=epsilon: epsilon_constraint(x, eps)}]

        # 初始解
        #x0 = np.random.uniform(bounds[0, 0], bounds[1, 0],bounds[2, 1])
        x0 = [15, 100, 0]
        bounds = np.array([[15, 30], [100, 130], [0, 30]])
        # 使用SciPy的minimize函数求解单目标优化问题
        result = minimize(f1, x0, method='SLSQP', bounds=bounds, constraints=cons)

        if result.success:
            solutions.append(result.x)

    return solutions


# 参数设置
bounds = np.array([[15, 20, 25, 30], [100, 110, 120, 130], [0, 10, 20, 30]])  # 变量x的取值范围
epsilons = np.linspace(0, 6, 20)  # ε值的范围和步长

# 运行ε约束法
solutions = epsilon_constraint_method(f1, f2, bounds, epsilons)
# X轴标签
x_labels = ['y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7']
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False   # 正常显示负号
# 绘制折线图
plt.figure(figsize=(10, 6))
# 输出结果
i=0
x11 = np.zeros(20)
x22 = np.zeros(20)
x33 = np.zeros(20)
for sol in solutions:
    #plt.plot(epsilons[i], sol[0], marker='o', label='x1')
    x11[i]= sol[0]
    x22[i] = sol[1]
    x33[i] = sol[2]
    i=i+1
    print(f'Solution: x = {sol}, f1(x) = {f1(sol)}, f2(x) = {f2(sol)} ')

plt.plot(epsilons, x11, marker='o', label='x1')
plt.plot(epsilons, x22, marker='o', label='x2')
plt.plot(epsilons, x33, marker='o', label='x3')
plt.xlabel('epsilons值')
plt.ylabel('最优工艺参数')
plt.title('不同epsilons对应的最优参数')
plt.legend()
plt.grid(True)
plt.show()